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Abstract 



Motivated by the present interest in the heavy quark effective theory, we use the 
spectator equation to treat the mesonic bound states of heavy quarks. The kernel 
we use is based on scalar confining and vector Coulomb potentials. Wave func- 
tions are treated to leading order and energies to order 1/rriQ in the heavy-light 
systems, and order I/nig in heavy-heavy systems. Our results are in reasonable 
agreement with experimental measurements. We estimate two of the parameters 
of the heavy quark effective theory, and propose further calculations that may be 
undertaken in the future. 
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Motivated by the present interest in the heavy quark effective theory, we use the 
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experimentai measurements. We estimate two of the parameters of the heavy quark 
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I. INTRODUCTION 

Recently, there has been great theoretical interest in hadrons containing b 
and c quarks. This has stemmed largely from the realization that, in the formal 
limit when the mass of one of the quarks in a hadron is taken to infinity, symme- 
tries above and beyond those usually associated with quantum chromodynamics 
(QCD) arise. This realization has led to the development of the heavy quark 
effective theory (HQET) [1] [2] [3]. In the framework of this effective theory, cor- 
rections to the formal limit can be systematically included. One very important 
phenomenological consequence of this has been a number of attempts to extract 
V c b from experimental data, with little model dependence in the result. 

Despite the power inherent in HQET, there is still much that this effective 
theory can not tell us about the properties of heavy hadrons. As an example, 
HQET allows us to infer the absolute normalization of some of the form fac- 
tors necessary for describing the decays of hadrons with beauty to those with 
charm. We also know how to include, in a systematic way, corrections to these 
normalizations due to the finite masses of the b and c quarks, as well as those 
due to perturbative QCD effects. We can even deduce bounds on the slopes of 



these form factors at a particular kinematic point. However, we know nothing 
about the exact dependence of these form factors on kinematic invariants. As 
a second example, HQET leads us to the conclusion that the spectra of B and 
D mesons should be very much alike, modulo 1/mj and \jm c effects. However, 
this effective theory tells us nothing about the details of the spectra, such as the 
exact ordering of states, or their masses. In essence, HQET provides a framework 
for systematically extracting symmetry relations and the corrections to the for- 
mal heavy-quark limit but can predict neither the spectra of the heavy mesons 
nor the approach to the heavy-quark limit. Until we know how to solve non- 
perturbative QCD, the details mentioned above, along with many others, are the 
realm of models: such models continue to play a crucial role in our understanding 
of QCD. 

A model that is quite successful in predicting the mesonic spectra is the rel- 
ativised constituent quark model of Godfrey and Isgur [4]. Indeed, it was this 
model and its applications to weak decays that originally suggested the existence 
of heavy-quark symmetries which in turn led to HQET. This model provides 
relativistic kinematic corrections to the standard nonrelativistic quark model us- 
ing a linear confining potential and a color Coulomb interaction. Meson spectra 
calculated with this model are remarkably close to experimental masses in all 
flavor sectors. However, since one of the objectives of heavy quark theory is the 
calculation of weak decay amplitudes and form factors, it is necessary to use a 
relativistically covariant model. 

A covariant extension to the Godfrey-Isgur model can be constructed using 
the spectator or Gross equation [5], which has been used with some success in 
models of the nucleon-nucleon interaction [6], as well as in quark models of mesons 
composed of equal mass quarks and antiquarks [7]. This equation can be related 
to the Bethe-Salpeter equation by placing one of the intermediate-state particles 
on the positive-energy mass-shell. This has the advantages that the prescribed 
constraint on the relative energy is manifestly covariant and that in the limit that 
the mass of one constituent goes to infinity (the static limit), the wave equation 
reduces to the Dirac equation for the light particle [8]. This is a property of the 
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full Bethe-Salpeter equation that is lost when the infinite sum of contributions 
to the kernel is truncated. Clearly, the properties of the spectator equation make 
it ideal for studying the properties of heavy mesons at finite mass. 

In this article we use the spectator equation to construct a constituent quark 
model of heavy mesons. In particular, we will use the spectator equation as 
a basis for construction and expansion of the heavy meson spectra and wave 
functions in I/tjiq, where mq is the heavy quark mass. This allows us to study 
the heavy meson spectra in the approach to the heavy quark symmetry limit. By 
choosing a reasonable set of model parameters we are able to obtain a respectable 
fit to the observed heavy meson masses and to predict the approximate masses 
of heavy mesons which have not yet been observed. 

This article is organized as follows. In the next section, we describe the model 
that we use for heavy mesons, including the derivation of a wave equation from 
the spectator equation. In Section III, three methods of obtaining solutions of 
the wave equation are described, while in Section IV we display our results. In 
Section V, we present some conclusions. 



II. THE MODEL 

A. Qq and qQ mesons 

The spectator equation is most easily understood in relation to the Bethe- 
Salpeter equation. The Bethe-Salpeter vertex function for two bound fermions is 
represented by Fig. 1 and can be written as 

T(p,P) = iJ ^V(p,k;P)S^\k 1 ,m 1 )S ( p 2) (k 2 ,m 2 )T(k,P), (1) 

where p = \{p\ — p 2 ), k = \{k\ — k 2 ), V is the Bethe-Salpeter kernel and 
Sp\ki,rrii) is the free Dirac propagator for particle i. The Dirac indices are 
suppressed for simplicity. 

The spectator vertex function can be obtained from the Bethe-Salpeter vertex 
function by placing one of the fermions on its positive-energy mass-shell. For our 
model the heavy quark (particle 2) is placed on shell while the light quark (particle 
1) remains off shell. This is achieved by a replacement of the propagator 

S { F ] (k 2 , m 2 ) - -2th- Fn m2 6 (k°-^- + E(k 2 ,m 2 )) A+ (2 \k 2 , m 2 ), (2) 
£j(K 2 ,m 2 ) \ Z J 

where 
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FIG. 1. Feynman diagrams representing the equation for the Bethe-Salpeter vertex 
function. 



A +(2) (fc 2 , m 2 ) = « (2) (k 2 , 4, m 2 )«( 2 )(k 2 , 4, m 2 ), (3) 

and replacing p, k and k\ by the corresponding quantities p, k and k\ with 
particle 2 on mass shell. The on-shell energy is given by -E"(p, m) = \/p 2 + m 2 . 
The spectator vertex function is then 

mP) = j ( 3°* 3 E( C 2 2 m2 ) V & k; ^)4 1} (fci^i)A +(2) (fc2, m 2 )T(k, P). (4) 
Defining the spectator wave function as 

ips 2 (P,P) = S ( p{pi,m l )u i2) {p 2 ,s 2 ,m 2 )Y(p,P), (5) 
the wave function satisfies the wave equation 

d) _1 

sy (Pi, m i)^s 2 (P,P) = 

where 

V S2 , S >(P, h P) = «( 2 )(p 2 , s 2 ,m 2 )V(p, k; iV 2 )(k 2 , 4 m 2 ). (7) 

This wave equation is covariant and can be easily boosted from frame to frame. 
It is generally easier to solve the wave equation in the bound-state rest frame 
where the angular expansions of the wave function and potential are defined. In 
the rest frame P = (W, 0), pi = — p 2 = p, ki = — k 2 = k, p° = W — E(p, m 2 ), 
p^ = E(p, m 2 ), ki = W — E(k, m 2 ), and k 2 = -E"(k, m 2 ) where W is the bound- 
state mass. The wave equation can be written as 

[jW°(W - E(p, m 2 )) - 7 (1) • P - mi] ^ S2 (p, W) = 
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where 



V S2<s .(p, k; W) = « (2) (-p, 8 2 , m 2 )V(p, k; W)u^(-k, s' 2 , m 2 ). (9) 

Since we wish to examine the approach to the limit m 2 — ► oo, it is useful to 
rewrite this equation in a noncovariant form by defining 



^^^Eik)^ (p ' w) (10) 



U, a> ,.(p, k; W) = /—H*— V , a> ,.(p, k; W). ™ 2 (11) 
12 w E(p, m 2 ) 2 V E(k, m 2 ) 



to give 



[ 7 (i) ( W _ E {p, m 2 )) - V 1 ) • p - mi] * S2 (p) 

/||E^(p,k;^(k). 



(12) 



It is necessary to assume some form for the kernel V in order to expand about 
the infinite mass limit. Here we assume that the kernel is of the simplest form 
which can be reduced to that used in ref. [4]. We choose the kernel to be 

V( P , k; P) = V S (Q 2 ) + 7 W • 7 (2) K(Q 2 ), (13) 

where 

Q 2 = (k - p) 2 - [E(k, m 2 ) - E(p, m 2 f . (14) 

V V (Q 2 ) is a vector potential which is a color Coulomb interaction and the confin- 
ing force is the result of the scalar potential V S (Q 2 ). This choice of interaction 
assumes that the Lorentz gauge is used in the color Coulomb interaction. 

Using the explicit form of the Dirac spinors in (12) and the Dirac 7-matrices 
to reduce particle 2 to the Pauli spin space, and defining a wave function which 
is an operator in the Dirac space of particle 1 and the Pauli space of particle 2, 
* = E s < X.'*.', (12) becomes 



(jW°(W - E(p, m 2 j) - 7 W • p - mi) *(p) 
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d 3 k f (E(p, m 2 ) + m 2 )(E(k, m 2 ) + m 2 ) 



(2tt) 3 V 4E(p,m 2 )E(k,m 2 ) 

1 ^.p^-k \ y ^ 

(E(p,m 2 ) + m 2 )(E(k,m 2 ) + m 2 ) ) >v *> 
+T V (E(p,m 2 ) + m 2 )(E(k,m 2 ) + m 2 )J " W J 

+ y . ( /> ->f + ) } * (k , (15 , 

\(E(p,m 2 ) + m 2 ) (E(k,m 2 ) + m 2 ) / 
Expanding eq. (15) to order l/m 2 , we find 

7 W°(W - m 2 - - 7 (1) 'P-mij *(p) = 

^ V s (q 2 )+7 (1) °K(q 2 ) 



(27T) 

+ 77^7 (1) • (<t (2) <t( 2 ) • k + <t^ ■ p^)K(q 2 )) *(k). (16) 
2m 2 J 

where q = k — p. 

Eq. (16) can be Fourier transformed to coordinate space, multiplied from the 
left by and then rearranged to give the wave equation 

#tf(r) = WV(r), (17) 

where the hermitian hamiltonian is H = Ho + Hi with 

H = «W • tV + /3 (1) m! + (]^V s (r) + V v {r) + m 2 , (18a) 




(r)}, (18b) 



where r is the unit vector in the radial direction. 

Eq. (18a) is the Dirac equation for particle 1 with scalar and vector potentials 
plus the mass of the heavy quark, particle 2. The solutions of the Dirac equation 
with such a potential have been extensively studied. The operators 

{j (1 > 2 ,#U,S< 2 >} (19) 
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TABLE I. Values of I and I for various values of k 



I I 


K\ < 


h - 


i 

2 


Jl + 


i 

2 


Kl > 


h + 


1 
2 


h - 


1 
2 



are a set of mutually commuting operators which commute with Ho, where 



L + S (i), S (i) = I S (i) = I^W 1 ), A'(!) = /^(SW - i) and S( 2 ) = . 



The eigenstates of Ho can then be labelled by the corresponding set of quantum 
numbers {n, j\, rrij 1 , K\, 82). The wave equation associated with Ho can then be 
written as 



#0* 



(o) 



(r) 



W (0) ■ * (0) • 

nreiji ^nsijirajj s 2 



(20) 



where 



* (0) • (r) 



G ntJ1 (r) m n ( 



(21) 



with 



j 1 m jl ) Y imi (tt)x Sl , 



(22) 



and Xsi an d Xs 2 are the Pauli spinors for particles 1 and 2, respectively. The 
eigenvalue k\ = ±(ji + ^) can be any nonzero integer. The values of I and I 
associated with various values of K\ are displayed in Table I. 

Note that the zeroth order invariant mass wi°) „• is determined by n, Ki and 
ji, or equivalently by n, j\, and £. The parity of the Qq bound state is given by 

p = (-iy+ i . 

The first term on the right hand side of (18b) is the kinetic energy of particle 
2. Both the first and second terms on the right hand side of (18b) commute with 
the set of operators given in (19). However, the third term does not commute 
with any of these operators, but instead commutes with 



\ J 1 Jz 



(23) 



where J = + S^ 2 ) and V is the parity operator. The eigenstates of the total 
hamiltonian H = Ho + Hi can then be labelled by the set of quantum numbers 
{n,J,Mj,P}. 
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The eigenstates and eigenenergies of the hamitonian H can be calculated 
directly. However, the objective of the calculations presented here is to produce 
wave functions which can be used in the calculation of form factors and decay 
constants as an expansion in powers of the inverse of the heavy quark mass 1112- 
In order to maintain consistency in this expansion, the masses and wave functions 
should be calculated perturbatively. The first order correction to the quark bound 
state mass is given by 

W$ P = j rf 3 ^,^^)^^,^^), (24) 

where 



J M j)^nl hmjl ,M (25) 



The bound state mass to first order is 



W n j P = wi°l n +W^ p . (26) 

The scalar and vector potentials in the calculations presented here have the 
form 

V s (r) = br + c, (27) 

A 3 

Vv(r) = --J2^f( 7t r). (28) 

8 = 1 

The vector potential is, as in ref. [4], based on a parametrization of the running 
QCD coupling constant. 



B. QQ mesons 

The situation for mesons made of a heavy quark and the corresponding an- 
tiquark is somewhat more complicated. The problem is that the prescription 
of placing particle 2 on mass shell in the Bethe-Salpeter vertex equation (1) to 
obtain the spectator vertex equation (4) is clearly asymmetrical. This results 
in a spectator vertex function which is no longer an eigenfunction of the charge 
conjugation operator. The solution of this problem is to construct a set of cou- 
pled equations for the vertex functions which have either particle 1 or particle 
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2 on mass shell [7]. These equations have been solved in ref. [7] for (/(/-systems 
containing only light quarks. 

However, since we are interested in expanding about the infinite mass limit, 
this additional complication is not necessary and a hamiltonian with leading 
1/rriQ corrections can be constructed from (4). The starting point is the spinor 
decomposition of the Dirac propagator of particle 1 in the meson rest frame 



5'^! ) (A;i,m (? ) 



m Q 



E{k,m Q ) ^7 



E 



wWfk, s[ , m Q )«( 1 )(k, s[ , m Q ) 



W -2E(k,m Q ) + ir) 



v (1) (-~k, s[, mq)v^\-k, s[, m Q ) 
W - IT] 

Using eqs. (29) and (3) in eq. (4), we can write [9] 

d 3 k 



m Q 



-V(p,k;P) 



where 



'(k) 



and 



) 3 E(k,m Q ) 

2 

( 1 )(k, S ' 1 ,m Q ) U ( 2 )(-k,4,mg)*g,(k) 

+v^(-k, s[, m Q )«( 2 )(-k, 4, m )*g, (k)) , 

m Q uW(k, s[, m Q )u( 2 H-k, 8 ' 2 ,m Q )T(k, P) 
E(k,m Q ) W-2E(k,m Q ) 

m Q vW(-k, s[, m Q )u( 2 H-k, .s' 2 ,m Q )T(k, P) 



E(k,m Q ) W 
Multiplying the terms of (30) to the left respectively by 

F/ m<? — T« (1) (P, s l ,m Q )u {2) {--p, s 2 , m Q ) 



and 



— — rf (1) (-p, si, m Q )M (2) (-p, s 2 ,m Q ), 
E(p,m Q ) 



(29) 



(30) 



(31) 



(32) 



(33) 



(34) 



means that eq. (12) can be rewritten as the pair of coupled integral equations 
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S 1> S 2 

+ ^:. a ;.' 1 ,.i(P' k ; wr )*ii( k )]> ( 35 ) 



and 

j3 



S 1> S 2 

+ U ^ S ^(P^ W )^ S ^)\> ( 36 ) 



where 



m Q 



K + s 2 ;s[,s' 2 (P> k 5 ^) ^ (pj mQ )E(k, m Q ) 

x«( 1 )(p, Sl )«( 2 )(-p, S2 )T/(p, k; W>W(k, S ;)«( 2 )(-k, 4), (37) 

U SltS2;s[A (p, k; = ^ mg)i?(k; mg) 

x«( 1 )(p, Sl , m Q )«( 2 )(-p, s 2 )T/(p, k; W>W(-k, 4)« (2) ("k, 4), (38) 



CA »i,» a ;»i,»^ P ' k; ^ = E(p,m Q )E(k,m Q ) 



xvW(-p, Sl )u^(-p, s 2 )V(p, k; W>W(k, S ;)« (2) (-k, 4), (39) 



U >i,>2;>' 1 ,>'S p ' k '' W ^ = E(p,m Q )E(k,m Q ) 



xv (1 H-p, Sl ,m Q )u (2 \-p, s 2 )V(p, k; W)v (1 \-k, 4)« (2) (-k, s' 2 ). (40) 

These coupled equations can then be reduced to the Pauli spin space and 
expanded in powers of 1/mg. In this case, only contributes to order 

l/rnq. Defining a wave function which is an operator in the spin spaces of both 
particles as 

* = £*«iX« a *S a ' (41) 
»i.» a 

eq. (35) becomes 

V - 2m Q - ^ *(p) = | 70^(P> k)*(k), (42) 
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where 



u(p, k) = y s ( q 2 ) + K(q 2 ) - [(^'(q 2 ) + K'(q 2 )) (k 2 - p 2 ) 2 

+ K (q 2 ) (p 2 + k 2 + <tW • p^ 1 ) • k + trW • p^ 2 ) • k) 
+ K (q 2 ) (p 2 + k 2 - <tW • po-^ • k - <t( 2 ) • p^ 2 ) • k) 

- K(q 2 ) (<t (1) <tW . k + ■ po-^) • (V 2 )<t( 2 ) • k + ■ p<r^] . (43) 

Eq. (42) can then be Fourier transformed to coordinate space to extract the 
hamiltonian 



with 



H — Hq + Hi, 



(44) 



where 



Hi — H c + ffhyp + H so + Hsr + ffyR, 



V 2 

H = + V s (r) + V v (r) + 2m Q , 

m Q 



(45) 



(46a) 



H c 



1 f 1 



m Q | 4 



[W s (r)] - [V v (r) - V s (r)]V 2 + [V^r) - V^r)} 



dr 



(46b) 



H 



hyP " ml 



1 [ 1 

2 



+ [V 2 K(r)] 



S -rS - r- rS : 



3^ 2 



(46c) 



ffs, 



2m^r 



[3!/'(r) - T/'(r)] S • L, 



(46d) 



H 



S(V)R 



Ami 



[V 2 ,[vVs (v) r(x)]] 



(46e) 
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and S = + S^ 2 ). Here _Fs(v)r( x ) is the Fourier transformation of 

^Ks(u)(q 2 )/^q 2 - F° r our choices of V s (r) and V v (r), we find 




2m qr* 



2 ^X) a ^ e_7<r 10 ^ " 4 ^ r ' + 8 ^ r 7> 7T- 4 -r>2 

offlgVi V V or r or or A 



(46g) 



Eq. (46a) is the nonrelativistic hamiltonian for equal mass quarks in scalar 
and vector potentials. H c contains central and orbital contributions, -ffhyp is the 
hyperfine interaction consisting of a tensor-force term and a spin-spin interaction. 
ff so is the spin-orbit interaction. ffsR and ffyR are scalar and vector retardation 
terms associated with the third term on the right-hand side of (43). Note that 
our spin-dependent interactions ffhyp and ff so have the same forms as those in 
many other quark models (see for example: [4,10,11]), but the spin-independent 
interactions do not. 

The spin-independent correction includes ff c , ffsR and ffvR- In these con- 
tributions, ffsR, ffvR and the term [V,'(r) — V^'(r)] in ff c are gauge depen- 
dent. ffsR and ffyR are from the second term in the expansion of V(Q 2 ) = 
^(q 2 ) — 5^2-^'(q 2 ) (k 2 — p 2 ) + 0{l/mq). Had we chosen the Coulomb gauge, 
these terms would not exist. Most other quark models do not include retarded 
interactions. (Ref. [12] gives another expression for the retardation effect.) We 
will show that with the scalar and vector potentials in (27) and (28), retardation 
contributions are comparable with the spin-dependent interactions. 

The operators {ffo, L 2 , S 2 , J 2 , 3 Z } where J = L + S, are a set of mutually 
commuting hermitian operators. The eigenstates of ffo can then be labelled by 
the corresponding set of quantum numbers {n, L, S, J, Mj}. The wave equation 
associated with ffo can then be written as 

tfo*i°i, JMj (r) = W<?MLjmA*), (47) 



where 



<°2sW*) = ^*,(n)> (48) 
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and 

y?sA Q )= J2 (LM l SMs\JMj)Y LMl (Q)\SMs) (49) 

M L , M s 

is the spin spherical harmonic. 

The hyperfine interaction (46c) mixes states with AL = ±2 for S = 1. As 
a result, L is no longer a good quantum number for solutions of the complete 
hamiltonian. However, these states have the same parity and charge quantum 
numbers since P = (— 1) L+1 and C = (— 1) L+S for . The first-order correction 
to the mass can then be written as 

W$ PC = J d^l\ JM ^)H^lsj M M 

= E C + E hyp + E so + E SR + E VR . (50) 
where P = (— \) L+1 and C = (— \) L+S . The bound state mass to first order is 

W nJPC = W$ + W$ PC (51) 

One may also include an annihilation term in the hamiltonian. However, this 

term first appears at order — f- [13] [4], while in our model the leading spin- 

dependent effects are of order ^-j-. Since a s is small in the heavy quark system 

(a s (m 2 ) ~ 0.35 and a s (m 2 ) ~ 0.22), we expect the annihilation effects on QQ 
spectra to be small. 



III. SOLUTION OF THE WAVE EQUATIONS 
A. Qq sector 

The Dirac equation (20) can be reduced by using the explicit forms of the 
zeroth order wave function (21) and the Dirac matrices a and [3 along with the 
identity 

< r( 1 )-f^(fi) = -^(fi) (52) 
to extract the coupled radial wave equations [14] 



dGni d J ; ir) + y G "th(r) = (mi + V s (r) - V v (r) + E^jF^r), (53) 
- ^F nln {r) = (mi + V s (r) + V v (r) - E^JG^r), (54) 
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where 



E ( °i = wi% n -m 2 . (55) 

We have obtained three separate numerical solutions of these coupled 
equations using two different techniques, direct integration and the matrix 
diagonalization- variational technique. 



1. Direct Integration 

This approach uses stepping techniques to obtain solutions to the differen- 
tial equations. Such techniques are much more efficient if any large asymptotic 
damping of the radial wave functions can be extracted and reduced radial wave 
equations can then be integrated. The scale of the asymptotic variation of the ra- 
dial wave functions is determined by the string tension b appearing in the scalar 
potential (27). Defining a dimensionless radial variable p = & 1/,2 r, and deter- 
mining the asymptotic behavior of the radial wave functions, the reduced wave 
functions g(p) and f(p) are defined in terms of G and F by 

G(r)=g(p)e-^ 2+ ^\ 

F(r) = f(p)e-i (p2+ ^ p \ (56) 

where j = 2(mi + c)/^ 1 / 2 , and c is the constant shift in the scalar potential. 
Coupled equations for the reduced wave functions that result are 

{t p ~ 9 ~ \ + 7") g(p) = {a+ + 9 - Vv(p)) (57) 

(J^-p-l-j)f(p) = (<*-+P+V v (p))g(p), (58) 

where V v {p) = K(r)/& 1/2 , a± = £ ± e and e = E^Jb 1 / 2 . 

In order to integrate the differential equations it is necessary to know the 
values of the functions and their derivatives at some point and then to have a 
stepping algorithm that predicts the values of the functions and their derivatives 
at subsequent points. The values of the functions and their first derivatives 
at p = are obtained by construction of a series solution for the functions for 
small p. An adaptive Runge-Kutte routine [15] is used to integrate the differential 
equations for increasing values of p. Energy eigenvalues can be found by adjusting 
the value of the energy until the functions have the correct asymptotic behavior 
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as determined by an asymptotic expansion of the functions at some large finite 
p. This process of finding the eigenenergies is called the shooting method [15]. 
In the calculations shown here, the accuracy of the eigenvalues is increased by 
integrating up from p = and down from some large finite p to some intermediate 
point where the values of g(p) and f(p) are required to match. 

A second variation on this method is to use the reduced radial wave equations 
(57) and (58) to eliminate f(p) to obtain a second order differential equation for 
g(p). This equation can then be integrated in a manner similar to the Shrodinger 
equation for the QQ sector. 



2. Variational Method 



The starting point for the 'variational' solution of eqs. (53, 54) is the pair of 
equations 



E^&tir) = (™i + V. + K)^(r) + -^—^(r) 



r J dr 



J nlj- 

= (V v - mi - V^r) + ^±%r) + 

^ (r) = ^W, Qltir) = (59) 

The functions T and Q are expanded in a set of orthonormal basis functions 
<f>\(r/ e ) 

N 



s&W = £«"#(»•/<?), 

8 = 1 

N 

^W = EA"«, (60) 



with 



drr 2 ^*(r/e)^(r/e) = 6 ltk . (61) 

o 

q is the size parameter of the wave functions, and is used as the variational 
parameter in this calculation. 

Substituting the expansion of eq. (60) into eq. (59), multiplying by 4>^*^(r/g) 
and integrating, leads to the set of equations 
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N , . 

.-1 \ / 1 



N 

> ' kl,il 

N , . 

= Y J (Vv(r)-m 1 - V s (r) ) ft 

i = l \ I kl,il 

' / kl,il i = i \ / kl,il 



8 = 1 

where we use the symbolic notation 



U(r)\ = ^ drr 2 4>1:(r/Q)Hr)4>U r /e)- (63) 

\ / kl 1 ,il 2 Jo 



The two sets of equations represented by eq. (62) can be combined into the single 
eigenvalue equation 



(m + V s (r) + V v (r)-E) 



d 



' (V v (r)- r m-V%)-E)){p)-°- (64) 

The size of the matrix in eq. (64) is 2N x 2N. Solutions to the eq. (59) 
are obtained by varying the wave function size parameter g, diagonalising the 
matrix in eq. (64) for each value of g, and searching for stationary points in the 
eigenvalues as functions of g. In principle, if the size of the expansion basis N is 
taken to oo, solutions obtained in this way would be exact and independent of 
g. In practice, the procedure outlined above is carried out for finite N , and N is 
increased until the eigenvalues are largely independent of g, for some reasonable 
range in g. With this method, the lower N eigenvalues obtained correspond to 
negative energy states, while the higher N eigenvalues are those of interest for 
this problem. 

For this problem we have used harmonic oscillator wave functions for the 
expansion, with N = 10 and N = 20. We compare the numerical solutions 
that we obtain using this procedure with those that are obtained using the other 
previously described methods. As expected, the variational solutions are better 
for N = 20, and the eigenvalues are within 1% of those obtained by solving the 
equations by the methods described in the previous subsection. 
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0.0 1 1 1 1 1 1 o.o I 1 1 1 1 1 

0.0 0.5 1.0 1.5 2.0 2.5 0.0 0.5 1.0 1.5 2.0 2.5 

l/p(GeV) l/p(GeV) 



FIG. 2. Energy eigenvalues as a function of ±, for N = 10 and N = 20 . 

B. QQ sector 

Using eq. (48) in eq. (47) and defining p = & 1/,2 r, the differential equation for 
the radial wave function is 

u„l(p) = eu nL {p), (65) 

where p = mq/b 1 ! 2 , e = {W^ - 2m Q - c)/b^ 2 and V v (p) = V v (r)/b^ 2 . 
Determining the asymptotic behavior of the radial wave function, the reduced 
radial wave function g(p) can be defined by 

u nL (p) = g(p)e-^(i^-<^ (66) 

The appearance of fractional powers of p in the argument of the exponential 
function in (66) leads to coefficients with fractional powers of p in the differential 
equation for g(p). This complicates the expansion of the reduced radial wave 
functions for small and large values of p. It is, therefore, convenient to define the 
variable £ = p 1 ! 2 . The differential equation for can then be written as 



I d 2 L(L + 1)\ — , . 
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+ AL(L + 1) + p-(sti + 2£ d ) - pe z e + 4^ 4 ^(r) 



g(0 = o (67) 



This equation can be used to develop expansions for small and large £ to provide 
boundary conditions for numerical integration of the differential equation. 

Since the Runge-Kutte method is designed to integrate systems of coupled 
first-order differential equations it is necessary to reexpress the differential equa- 
tion (67) as the coupled pair 



■^9(0 = f(0, 



(68) 



and 



-e^ + (t-2sp H 2 - 



Ape 



f(0 



+ AL(L + 1) + p-(sti + 2£ d ) - pe 2 e + AptY^) 9(0 = 



(69) 



This system can then be solved by Runge-Kutte integration and intermediate- 
point shooting techniques. 



IV. RESULTS 



Once the zeroth-order solutions are found, the perturbed energies can be 
calculated using (24) and (50). The masses associated with the bound states are 
given by (26) and (51). These depend on the quark masses m u , m s , m c and mj 
as applicable for each meson; the parameters of the scalar potential (27) b and 
c; and the parameters of the vector potential (28) a; and ji for i = 1,2,3. The 
model contains a total of twelve parameters. In obtaining the results shown here, 
the vector potential parameters 

a-j = 0.15, «3 = 0.2, 

7i = 0.5, 7 2 = 1.581, 7 3 = 15.81, (70) 

are fixed at the same values as given in ref. [4]. The remaining vector potential 
parameter u\ is reexpressed as 

ai = a crit - a 2 - a 3 . (71) 
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TABLE II. Parameters ol the model. 



parameter value comments 

a crl t 0.674 limiting value ol a s 

b 0.180 GeV 2 string tension 

c 0.02 GeV see eq. (27) 

m u 0.258 GeV 

m 3 0.400 GeV 

m c 1.53 GeV 

m b 4.87 GeV 





TABLE III. 


Fitted meson spectra for Qq mesons. 








Mass (GeV) 




Meson 


J p 


theory 


experiment 


D 


0" 


1.85 


1.87 


D* 


1" 


2.02 


2.01 


Di 


1 + 


2.41 


2.42 


Dt 


2+ 


2.46 


2.46 


B 


0" 


5.28 


5.28 


B* 


1" 


5.33 


5.33 


D s 


0" 


1.94 


1.97 


D* s 


1" 


2.13 


2.11 


B s 


0" 


5.37 


5.38 


B* s 


1" 


5.43 


5.43 



Experimental values are quoted [16] to the nearest 10 MeV due to ambiguities in 
assigning the calculated values to specific charge states. 



where a cr a is the value of the running coupling constant at Q 2 = as 
parametrized in ref. [4]. 

u cr it and the remaining model parameters are adjusted to fit the masses of a 
selection of mesons. The resulting values are listed in Table II. The fitted meson 
spectra for the Qq sector are listed in Table III and the fitted meson spectra for 
the QQ are listed in Table IV. Additional states which were not used in the 
fitting procedure were calculated and a detailed discussion of the results for the 
Qq and QQ is presented in the following two subsections. 
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TABLE IV. Fitted meson spectra for QQ mesons. 



Meson 


T VC 


theory 


Mass (GeV) 

experiment 


Vc 


T 

_+ 


3.00 


2.98 


J/f/>(lS) 


1 — 


3.10 


3.10 


XcO 


0+ + 


3.44 


3.42 


Xcl 


1++ 


3.50 


3.51 


Xc2 


2++ 


3.54 


3.56 


J/i>(2S) 


1 — 


3.73 


3.69 


T(15) 


1 — 


9.46 


9.46 


Xbo(lP) 


0+ + 


9.85 


9.86 


Ybl(l-P) 


1++ 


9.87 


9.89 


Xta(lP) 


2+ + 


9.89 


9.92 


T(25) 


1 — 


10.02 


10.02 


Xto(2P) 


0+ + 


10.24 


10.24 


Xti(2P) 


1 + + 


10.26 


10.26 


Xb2(2P) 


2+ + 


10.28 


10.27 


T(35) 


1 — 


10.39 


10.36 



A. Qq sector 

For the Qq sector, the zeroth-order eigenenergy E^^ = W^°J i j i — is 
independent of the heavy quark mass, as would be expected in the heavy quark 
limit, where the heavy quark should act as a static source. The zeroth-order 
spectrum depends only on the light quark mass. The first-order correction to the 
mass wj^jp is proportional to l/m? and splits each of the unperturbed states. 
These features are illustrated in Fig. 3 which shows W^°J i j i — for a u quark as 

solid lines and W n jp — rri2 = H^kIji + ^ijp — m 2 with a c quark as the heavy 
quark (dotdashed lines) and with a b quark as the heavy quark (dashed lines). 
Fig. 4 is a similar spectrum where the light quark is now an s quark. 

Note that to zeroth order the ordering of the j\ = I ± 1/2 states is reversed 
for the I = 2 states in comparison to the 1=1 states. This phenomenon, called 
multiplet inversion, has been predicted [17] for Qq mesons with rri2 ^> m\. It 
results from the dominance of the Thomas-precession over the spin-dependent 
forces in this limit. 

For the states presented here, the root mean square momentum of the zeroth- 
order wave function is approximately 0.9 GeV. Clearly, both u and s quarks are 
very relativistic. In addition, it is possible to obtain some sense of the convergence 
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1.4 
1.3 
1.2 
1.1 
1.0 
0.9 

> 

« 0.8 

O 

0.7 
0.6 
0.5 
0.4 
0.3 
0.2 



W ,0, -m 2 for mesons with u quark 
W cu - m 2 
W h „ - m. 



l=0j=l/2 l=lj=l/2 



l=lj=3/2 l=2j=3/2 



l=2j=5/2 



FIG. 3. This figure shows W — m,2 for bu and cu to the zeroth order and to the first 
order, h and ji are the quantum numbers for orbital angular momentum and total 
angular momentum of the u quark. The states have been labelled as J F . 
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1.4 
1.3 
1.2 
1.1 
1.0 
0.9 
4J0.8 

o 

0.7 
0.6 
0.5 
0.4 
0.3 
0.2 



W (0| -m 2 for mesons with s quark 
W C! - m 2 
W h , - m. 



l=0j=l/2 



l=lj=l/2 



l=lj=3/2 



l=2j=3/2 



=2j=5/2 



FIG. 4. This figure shows W — m,2 for bs and cs to the zeroth order and to the first 
order, h and ji are the quantum numbers for orbital angular momentum and total 
angular momentum of the u quark. The states have been labelled as J F . 
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3.60 



3.20 



2.80 



2.40 



3.37 



3.13 



2.98 3.07 



2.50 



2.71 



2.62 



3.16 



3.19 



2.76 2.78 



2.74 



3.24 



2.78 



3.30 



2.89 



3.38 



3 24 3.26 ^=^= 3 34 

3.17 _ =^^= 3.20 3.29 



2.27 



2.41 
2.40 



2.46 



3.41 
3.39 



3.46 
3.41 



3.00 3.03 ^03_ 
3.01 



2.94 



2.00 



2.02 



1.85 



data 

this work 



1.60 1 1 1 1 1 ' ' — ' t 

0" 1 2" 3" 4 1 2 3 T 4 T 

FIG. 5. cu spectrum. In this figure, solid lines represent the results of our calculation 

for the masses of cu mesons, W, to the first order in the perturbation; dotted lines 

represent the data. 



of the p/m expansion for the corrections to the infinite-heavy-quark-mass limit 
since Prms ~ ^ while Prms ~ ^ Therefore, the higher-order correction that 
are neglected here should be considerably larger for the the c quark than the 
b quark. Indeed, this problem will become worse with increasing n since p rms 
should increase with increasing n. This is seen in the shift of the 0~ states relative 
to the unperturbed states which increases with n. 

Figs. 5 to 9 show predictions for the masses of Qq mesons, W, to first order 
in the perturbation (solid lines). In the spectra for mesons with u and s quarks, 
the available data are plotted for comparison as dotted lines. Ref. [16] has also 
listed states Dj(2.440) and D s j(2.573) with uncertain quantum numbers. We 
believe they correspond to the state 1 + (2.41) in Fig. 5 and the state 2 + (2.58) 
in Fig. 6 respectively. For the be mesons, calculated masses from [4] are plotted 
because no data exist at present. For the be mesons, p ™ s ~ 1. This shows that 
although the mass of the c quark is relatively large it is quite relativistic in this 
case. 
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3.60 353 ^58_ 

3.50 J ' 5J 



3.25 



3.36 ^jg_ = 3.46 

3 - M ^37~^- 3 - 41 

3.31 



3.42 3.52 



3.54 



3 20 3.28 — 3 . 16 

r 3.19 3.12 



3.09 3.01 



2.88 2.90 2.90 3.00 



2.82 



2.80 2 - 86 



2.73 



2.61 2.52 

2.40 

2.13 

2.00 

1.94 

1.60 



3.06 



2.58 



3.13 



data 

this work 



0" 1" 2" 3" 4" 1 2 3 4 

FIG. 6. cs spectrum. See caption of Fig. 5. 
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7.00 



6-60 6 . 52 



6.32 



6.38 



6.49 650 



6.55 



6.48 6-49 



6.20 



5.80 



6.21 



624 6-32 
5.97 5 - 98 5.97 



6.19 6.20 618 



6.36 



^.-^^12- 



6.18 



-^—-JW 



5.96 



5.65 



5.69 



5.71 



5-40 533 

5.28 — 



5.00 



data 

this work 



1" 2" 3" 4" 1 2 3 4 

FIG. 7. bu spectrum. See caption of Fig. 5. 
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7.00 



6.60 



6.20 



5.80 



5.40 



5.00 



, , . 6.65 

6 - 62 6.60 6.61 

6.49 6.49 6.50 6 53 - . - 6.59 6-60 

6.42 6.52 6.54 



6.47 



6.31 634 6.41 6M 6 . 30 6.31 6.29 



5.93 — — 
5.97 



5.43 
5.37 



fi .07 6.08 6.08 
6.07 



6.17 



5.75 



6.21 



6.20 6.22 



5.80 5.82 
5.79 



6.28 



data 

this work 



1 2" 3" 4" 1 2 3 4 

FIG. 8. bs spectrum. See caption of Fig. 5. 
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8.00 



7.60 7.55 



7.20 



6.85 



6.80 



7.64 
7.58 
7.35 
7.28 



6.90 



7.65 


7.70 


7.71 


7.71 




7.70 






7.66 






7.43 


7.44 


7.43 - 


7.36 




7.43 






7.37 




7.03 


7.04 




7.10 - 



7.02 





7.80 


7.81 


7.85 


7.74 






7.80 


7.81 




7.76 




7.54 


7.55 


7.60 






7.47 




7.54 
7.25 


7.55 


7.46 


7.48 
7.24 




7.15 




7.24 


7.25 - 




7.16 







6.68 



7.14 



6.76 



6.73 



6.40 



6.26 



6.00 



6.34 



Godfrey and Isgur 
this work 



1 2" 3" 4 1 2 3 4 

FIG. 9. be spectrum. See caption of Fig. 5. 
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In these figures, the results are in good agreement with the data, which vindi- 
cates our choices of potentials and parameters. However, the calculated hyperfine 
splittings are all larger than in the data. The agreement is much better in the 
&-flavored mesons than in the c-flavored mesons. There are three possible reasons 
for this discrepancy. First, as has been mentioned earlier, this model is expected 
to work better for &-flavored mesons than for c-flavored mesons due to the more 
rapid convergence of the nonrelativistic expansion applied to the heavy quark. 
Secondly, these calculations do not include any effects associated with possible 
strong decay of the heavy mesons. The coupling to these strong decay channels 
will result in shifts in the meson masses as well as decay widths for heavy mesons 
above decay thresholds. These shifts will be greatest near the decay thresholds. 

The third possible reason for the large hyperfine splittings may have its origin 
in the parametrization of a s (r), particularly at small r. While many functional 
forms may be used for this parametrization, each form may be expected to lead 
to quite different 1/rriQ contributions, especially in the hyperfine term. This 
question is currently under investigation. 

The third term on the right hand side of (18b) has off-diagonal matrix ele- 
ments between states with j\ differing by unity and with I differing by either 
or 2. These mixings do not affect the spectrum to order but should result 
in shifts in some states at higher order in all of these systems. This should be 
particularly apparent for the 1 + states which are nearly degenerate to order 
for all Qq mesons calculated here. 

One very interesting aspect of this calculation is the mapping of our model 
onto the heavy quark effective theory, with a view to evaluating some of the 
parameters and dynamical quantities (such as universal form factors) of the ef- 
fective theory. While we do not endeavor to perform such a calculation for all 
such quantities here, some comments are merited. 

Although we have included all of the 1/rriQ terms that arise from the spectator 
equation, it is not clear that these correspond to all of the 1/rriQ terms of HQET. 
In particular, in the spectator equation, the heavy quark is treated as being 
exactly on its mass shell. In contrast, in HQET, the heavy quark is allowed 
to be slightly off its mass shell (via the equation = mqv^ + k^), and this 
leads to terms that may be absent from the formulation presented here. The full 
ramifications of this are also under investigation. 

Until this question is resolved, we dare not examine quantities that are inti- 
mately bound up in the 1/rriQ structure of the effective theory or the model. We 
can, however, examine quantities that depend only on the leading-order structure 
of the model, as we believe that this is a reasonably accurate representation of 
the effective theory. In particular, in the effective theory, one expects that the 
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heavy quark should act as a static color source. This very important feature is 
reproduced in the model, as the leading dynamical behavior is described in terms 
of a Dirac equation for the light quark. 

Two quantities of interest in HQET are A and Ai, which are defined by 

M M = m Q + A + O [— ) , 

\ m Qj 

(M(v)\h Q (iD) 2 h Q \M(v)) =2M M Ai. 

A is crucial for the effective theory, as it appears as the coefficient in the 1/rriQ 
expansion: the expansion coefficient is written as A/tjiq. A is, in essence, the 
contribution to the mass of the meson from the mass and kinetic energy of the 
"brown muck" . The left hand side of the second expression above is proportional 
to the kinetic energy of the heavy quark. The meson states in the bra and ket 
above are the leading order representation, and so correspond to our zeroth- 
order calculation. From our model, we obtain A = 0.45 GeV for the ground 
state pseudoscalar/ vector doublet, and Ai = 0.67 GeV 2 . These values are in 
reasonable agreement with other values in the literature [3]. Further aspects of 
the relationship of our model to HQET are discussed in the conclusions. 



B. QQ sector 

Figs. 10 and 11 show the spectra for cc and bb mesons as calculated with 
eqs. (44)-(51). As before, the calculated masses are shown as solid lines and the 
experimental masses as dotted lines. The DD and BB thresholds are shown as 
horizontal dotdashed lines across the Figs. 10 and 11 respectively. Ref. [16] has 
also listed states h c (lP) with mass 3.526 GeV and r) c (2S) with mass 3.590 GeV. 
We believe they correspond to the states 2 1 5'o(3.67) and l 1 Pi(3.51) in Fig. 10 
respectively. 

The bb spectrum is in quite good agreement with the data for the states lying 
below the BB threshold. The agreement deteriorates as the masses approach and 
cross the BB threshold. As argued in the previous section, this may be the result 
of the absence of coupling to strong decay channels. The agreement for the cc is 
less satisfactory. This may be an indication of the inadequacy of the truncation of 
the nonrelativistic expansion at order ■ In both cases the hyperfine splitting 

of the spin triplet states is too large. 

Since the hyperfine tensor interaction has non-zero off diagonal matrix ele- 
ments for states with spin 1 and with L differing by or 2, there should be 
mixings of states such as 3 Si with 3 D\ and 3 P'j with 3 i<2- These mixings do 
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4.80 r 

3 3 D,<4.59) 



4.40 



2^(4.22) 2'P 2 (4.24) 2 3 D;(4.24) 2D,(4.24) 

3 j S,(4.13) 3'Si(4J8) 

1 3 F;(4.06) l'F 3 (4.06) 1 3 F 3 (4.06) 1 3 F 4 (4.06) 



2 , P 1 <3.99) 2 3 P,(3.99) — 

4.UU - 2 3 P,(3.94) 2 3 P,(4.02) 

^ ,•»„.«., DD threshold ^ ^ ^ 

*J 2 3 S,(3.73) 

U^jfll L2'S„(3.67) 



3.60 

1 P,(3.51) l'PiO^O) 



3.20 

rSjO.10) 

l's (1 (3.00) 



2-80 



this work 



0" T 1" r 1" 2" 2" T 2" 3" 3 T " 3" 4 T 

FIG. 10. cc spectra. See caption of Fig. 5. 



not affect the spectrum to order but should result in shifts in some states at 

higher order in both the bb and cc spectra. 

Table V shows the individual contributions to the masses W of a number of 
bb states from W^°\ E c , -Ehyp, E so , Esr and EVr. The retardation contributions 
-Esr and £Vr are clearly gauge dependent since they would not appear in the 
Coulomb gauge. E c is also gauge dependent. These contributions may also be 
sensitive to the choice of quasipotential prescription. To this order Ehyp, E so 
should be independent of these factors. Note that the scalar and vector retarda- 
tion contributions are of opposite sign and therefore tend to cancel. However the 
sum of these contributions is comparable with -Ehyp and E so . The assumption 
that the scalar retardation potential depends only on the square of the exchanged 
four-momentum Q 2 is uncontrolled and it is possible to propose forms for this 
retardation potential which would eliminate the scalar term altogether. 
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11.20 



10.80 



> 10.40 

o 



10.00 



9.60 



9.20 



-6 3 S,(11.15) 



-5%( 10.93) 



3 3 D,(10.73) 



^ 3 S,(10.68) 



BB threshold 



2 3 D,<10.46) 2'P 2 (10.47) 2 3 D ; ( 10.47) 2 3 D 3 ( 10.47 ) 
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1^(9.87) l 3 P2<9-89) 



1 3 S,(9.46) 



data 

this work 



0" T 1" l T 1" 2" 2" T 2" 3" 3 T " 3 T 

FIG. 11. bb spectra. See caption of Fig. 5. 
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TABLE V. Zeroth order and various first order interaction energies in the bb spectrum 













(GeV) 








State 


w 


E c 






E SR 


E VR 


E SR + E VR 


l 1 So 


9.41 


9.5315 


-0.0602 


-0.0367 


0.0000 


0.0072 


-0.0297 


-0.0224 




9.46 


9.5315 


-0.0602 


0.0122 


0.0000 


0.0072 


-0.0297 


-0.0224 



2 1 S 
2 3 S 1 

3 1 S 
3 3 S t 

4 1 S 
4 3 S 1 

5 3 S t 

e 1 s 

6 3 Si 



10.00 
10.02 



10.37 
10.39 



10.66 
10.68 



10.91 
10.93 



11.14 
11.15 



10.0892 
10.0892 



10.4511 
10.4511 



10.7411 
10.7411 



10.9928 
10.9928 



11.2202 
11.2202 



-0.0708 
-0.0708 



-0.0839 
-0.0839 



-0.0992 
-0.0992 



-0.1162 
-0.1162 



-0.1345 
-0.1345 



-0.0192 
0.0064 



-0.0146 
0.0049 



-0.0125 
0.0042 



-0.0113 
0.0038 



-0.0105 
0.0035 



0.0000 
0.0000 



0.0000 
0.0000 



0.0000 
0.0000 



0.0000 
0.0000 



0.0000 
0.0000 



0.0175 
0.0175 



0.0302 
0.0302 



0.0447 
0.0447 



0.0608 
0.0608 



0.0781 
0.0781 



-0.0192 
-0.0192 



-0.0160 
-0.0160 



-0.0144 
-0.0144 



-0.0135 
-0.0135 



-0.0128 
-0.0128 



-0.0017 
-0.0017 



0.0142 
0.0142 



0.0303 
0.0303 
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V. CONCLUSION AND OUTLOOK 



We have constructed this model for heavy mesons based on a relativistic 
bound state equation, namely the spectator equation. The calculated spectra are 
in quite good agreement with the experimental data. The parameter values we 
have are reasonable, and comparable to other models of similar type. The model 
is derived by expanding the spectator equation in 1/Mq, where Mq is the mass of 
the heavy quark. This treatment is expected to work better for &-flavored mesons 
than for c-flavored mesons since in c-flavored mesons, v ~ -|c, but in &-flavored 
mesons, v ~ |c, and our results confirm this expectation. 

The retardation contribution to the QQ mesons, which is missing in other 
quark models, has a noticeable effect. Annihilation effects have been neglected, as 
they are suppressed by additional powers of a s (Mg), which is a small parameter. 

In addition to the questions currently being investigated (parametrization 
of a s (r), 1/rriQ terms), this work opens up many avenues of investigation. Of 
primary importance is the application of the model to decay processes of heavy 
mesons. In particular, the calculation of the Isgur-Wise functions that describe 
the semileptonic decays, not only for decays to pseudoscalars and vectors, but also 
to excited states, are of great interest. In HQET, these form factors are essentially 
the overlaps of the appropriately boosted wave functions. It will be interesting 
to see if this relationship between the form factors and the wave functions arises 
in the present model, and if so, how. In addition, the slope of the Isgur-Wise 
function for the elastic decays may also be calculated, and various HQET sum 
rules checked. 

The strong and electromagnetic decays may also be treated with the wave 
functions that we have. These are particularly interesting for the D* and D* 
states, as the former lie so close to the Dtt threshold, while the latter lie below 
the DK threshold, and thus decay radiatively. In addition, quantities such as 
meson decay constants may also be evaluated. 
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